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SUMMARY 

The paper presents a numerical method for the study of steady, transonic, 
turbulent, viscous flow through plane turbine cascades. The governing equa- 
tions are written in Favre-averaged form and closed with a first order model. 
The turbulent quantities are expressed according to a two-equation K-e model 
where low Reynolds number and compressibility effects are included. The solu- 
tion is obtained by using a pseudo-unsteady method with improved perturbation 
propagation properties. The equations are discretized in space by using a 
finite volume formulation. An explicit multistage dissipative Runge-Kutta 
algorithm is then used to advance the flow equations in the pseudo-time. 

First results of calculations compare fairly well with experimental data. 


INTRODUCTION 

As the efficiency levels of turbine engines continue to increase, the 
accurate prediction of blade performances becomes increasingly critical in the 
development and design process. Although numerical methods to solve steady, 
transonic, turbulent, viscous flows have been developed, efforts to apply 
these methods to the calculation of performance of turbine blades have so far- 
proved somewhat unsatisfactory. This is mainly due to the failure of present 
mathematical models to consistently simulate the complex phenomena inherent in 
these flows. Futhermore, additional problems arise due to the presence of 
numerical viscosity in the solution algorithm, sometimes of the same order of 
magnitude of the physical one. 

Previous calculations have been performed by using the simplifying assump- 
tion of local equilibrium conditions and then evaluating the turbulent viscos- 
ity coefficient via a generalized mixing length formula (ref. 1). These 
calculations have proved to lead to quite accurate predictions of blade pres- 
sure distributions at design point conditions, but in an analysis of blade per- 
formances, while qualitative behavior of loss generation has been correctly 
predicted, predictions of quantitative behavior must still be further improved. 
Furthermore, in performing calculations of separated flows, only qualitative 
results can be obtained by retaining the assumption of local equilibrium condi- 
tions, as clearly represented by backward facing step flow calculations. 


*Uni vers i ty vi si tor . 



The work presented here is focused on developing low-Reynolds number and 
compressible K-e turbulence models for solving the problem of correctly pre- 
dicting blade performances. The method includes low Reynolds number terms, so 
that the equations are valid all over the laminar, transition, and turbulent 
region (ref. 2). Furthermore, the method includes a density gradient term to 
better simulate variable density effects. 

The solution is obtained by using a pseudo-unsteady method with improved 
perturbation propagation properties (ref. 3). The equations are discretized 
in space by using a finite volume formulation. An explicit multi-stage dissi- 
pative Runge-Kutta algorithm is then used to advance the flow equations in the 
pseudo-time. Multi-stage schemes for the numerical solution of ordinary dif- 
ferential equations are usually designed to give a high order of accuracy, but 
in a pseudo-unsteady solution these schemes are selected only for their proper- 
ties of stability and damping. The enhanced stability of a four-stage scheme 
allows one to considerably reduce the numerical viscosity of the dissipative 
terms . 
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SYMBOL LIST 

constant 
Courant number 
diffusion vector 
total specific energy 
specific internal energy 
flux tensor 
unknown vector 
function 
total enthalpy 
optimization matrix 
i denti ty matrix 

index of the spatial discretization 
turbulence kinetic energy 
index of the multistage algorithm 
length scale of turbulent motions 
Mach number 
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m index of the temporal discretization 
N unit outward normal 

NV updating rate 

P production of turbulence kinetic energy 

Pr Prandtl number 

p pressure 

Q work due to turbulence 

q heat flux vector 

R low Reynolds number term 

Re Reynolds number 

S source vector 

Sc Schmidt number 

T residual 

t time 

tu turbulence intensity 

V velocity 

v volume 

x axial coordinate 

Y tangential coordinate 

a flow angle 

3 perturbation speed 

y specific heat ratio 

Sr characteristic volume dimension 

St time step 

S£ surface area 

e turbulence kinetic energy dissipation rate 

C kinetic energy loss coefficient 
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9 multistage scheme coefficient 

k thermal conductivity 

y viscosity coefficient 
p density 

E boundary of the fixed volume 

x viscous stress tensor 

Q artificial viscosity coefficient 

Subscripts 
i inviscid 

is isentropic 

1 laminar 

t turbulent 

v viscous 

0 total 

1 inlet 

2 outlet 


GOVERNING EQUATIONS 

The unknown vector f is the solution of a system of conservation equa- 
tions. This system is written Favre-averaged , dimensionless, vector, integral 
form as fol lows 


** p p 

N . <F. + F )dl = 

J J IV t 

E 

The basic dependent variables are density, 
conservation equations read as follows 


J J s dv 

v 

velocity, and energy. 


Their 



fp * V \ 

F. = p • V V + p • I 

Vp • E • v + p • I • vy 
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F v = r x 

V-t • V + q J 

s-(S ) 

Vp * e - P + Qy 

where from the equation of state of a perfect gas 

P = ( y - 1 ) • p • e 

wi th 


2 2 

e = (e - V /2 ) = (h - V /2)/y 

The stress tensor x is the sum of a laminar and a turbulent part, where 
the latter is expressed according to a classical eddy viscosity concept, i.e. 


x = -2/3 • (p/Re • div(V) + p • K) • I + 2 • p/Re • (def ( V) ) 


where 


P = + P t 

Similarly, the heat flux vector is given by 

q = ~Y • k • grad(e)/(Re • Pr) 


k = + kt 

The turbulent viscosity coefficient is expressed according with the 
Prandtl-Kolmogorov formulation 

2 

p. = C • f • p • Re • K /e 

t p p 

while 

k » p 

The turbulence variables are the turbulence kinetic energy and its dissi- 
pation rate. Their conservation equations are written in the following low 
Reynolds number and compressible form 


f = 
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(? • V . 
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P-Q-p*e + R^ 

(P - Q) • C , • j , • c/K - C - • f , • p • e 2 /K + R 
el J el e2 J e2 K e 

where the production of turbulent energy from the mean flow energy P and the 
work due to turbulence Q are given by 

P = p t - (def(V) : def (V) )/Re 



0 = Cp • ^t * (grad(p) • grad(p)/p 2 )/Re 

the low Reynolds number functions are given by (ref. 2) 

f = exp(-3 .4/(1 . + Re. /50.)> 
p t 


f 


el 


1 .0 


f e2 = 1.0 - 0.33 • exp(-Re 2 ) 

R e = -2. • p a • (grade(c)) 2 /Re 
R k = -2. • p • (grad(K) ) 2 /Re 

where 

Re^ = p • K 2 * Re/(p Jt • e) 
and the model constants are assumed as 

C , = 1 .43 C , = 1 .92 Sc =1.3 Sc. = 1 .0 

Cl C C- G l\ 

C = 0.09 C =1.0 

M P 

The above two equation model does not take into account the preferential 
damping of velocity fluctuations in the direction normal to the wall, but it 
is quite general and it is useful in laminar, transition, and turbulent 
regions. Furthermore, the model adopts the assumptions and approximations 
which are normally used for constant density flows, by retaining the gradient 
diffusion model to be rewritten in the density weighted form without any 
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explicit account being taken of density fluctuations. However, the introduc- 
tion of the compressibility term Q allows a partial consideration of varia- 
ble density effects. 

Along the inflow boundary, the total pressure, total density, flow angle, 
inlet turbulence level, and length scale of the turbulent motions are speci- 
fied, while the Mach number is extrapolated from the interior. Along the out- 
flow boundary, the static pressure is specified, and all the other variables 
are extrapolated. Along the solid boundaries, the no slip condition requires 
the vanishing of velocity, turbulence kinetic energy, and turbulence kinetic 
energy dissipation rate, the latter intended to be the modified quantity used 
in the conservation equations. Furthermore, for adiabatic flows, the specific 
internal energy gradient normal to the wall is set equal to zero. The density 
is finally extrapolated by the interior. 


NUMERICAL SOLUTION 

The proposed equations are solved in two-dimensional geometries by using 
a pseudo unsteady method with a finite volume, dissipative, explicit discreti- 
zation. The solution of the steady equations is obtained as the asymptotic 
solution of the following artificial unsteady equations 
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(If) • H" ] dv + 

\9t/ p J 


N • (F. + F )dE = 

1 V t 




These unsteady equations are generally constructed in order to obtain the bet- 
ter convergence rate, obviously providing that the steady state solution is 
not altered. 

From the identity between the convergence process and the elimination 
process of the initial perturbations to the steady solution, the convergence 
parameters are determined in order to improve the perturbation propagation or 
damping. We use (ref. 3) 


H = 6r/(3 
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f 2 = max(f 1 + 1 . , C ) 


wi th 


- 1/2 


M = V • (y • P/p) 

Cp is a small positive number, and 

^max " m1n (li r V 


(i ] = V • {(1. + M 2 )/2. + [((1. - M 2 )/2.) 2 + l.] ,/2 } 


H 2 = V 


(1 . + 1 . /M) 


These convergence parameters produce an improved ratio for subsonic flows 
between the speeds of the fastest and slowest perturbation (ref. 3), and there- 
fore result in an improved propagation. 

These equations are then discretized in space by using a finite volume 
discretization. The mesh is nonorthogonal and curvilinear, conforming to the 
boundaries of the domain, with lines intersecting at arbitrary angles, prop- 
erly refined where high gradients are expected to occur. The discretization 
nodes, located at the intersection of these lines, are the centers of hexago- 
nal control volumes, obtained by connecting the six surrounding nodes. A 
sample computational domain and the hexagonal control volume are shown in 
figure 1 . 

The discretized equations are written as follows (refs. 2 and 4) 


8f H • 
3t " P 


- 

' 6 

£ (F ,,j ♦ F v,i> • N J ' 

/♦*( 

, 

Lj-1 J 



where the subscript j refers to every face of the finite volume, 
tization is second order accurate on smoothly varied meshes. 


The discre- 


te equations are finally discretized in time by using an explicit, dis- 
sipative discretization. Let the previous equation be rewritten with the add i 
tion of a dissipative term as follows 


3f 

at 


= T(f ) + D(f ) 


where T represents the residual and D is the dissipative term. An explic- 
it k-stage Runge-Kutta algorithm, based on the work of Jameson (ref. 6), may 
be written as follows 


f(0) = f m 


f (l) = f (0) - Qj • St • [T(f (0) ) + D(f (0) ) ] 
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• St • [ T ( f 


f (k-l) = f (0) _ 0 


k-1 


<k_2) ) + D(f (k 2) )] 


: (k) = f (0) - St • [T(F (k n ) + D(f (k 1} )] 


^m+1 _ ^(k) 

A four-stages scheme, with the standard coefficients 

©! = 1/4 0 2 = 1/3 0 3 = 1/2 

has a Courant limit CFL = 2.8 (refs. 6 and 7). If the dissipative part is 
evaluated only once, the scheme has a slightly reduced Courant limit CFL = 2.6 
(ref. 6), but it is computationally more efficient. 

In order to further improve the computational efficiency, the dissipation 
is corrected and the viscous terms are evaluated only at fixed iterations. 

Then the k-stage scheme can be conveniently rewritten as follows 


f (1) = f (0) - e. 


f (k -i } = f (o> _ 


k-l 


f (k) = f (0> 


^(0) _ ^m 

• St • [T. (f (0) ) + Ty(f m *) + D(f <0) ) - Q • D(f m * ) ] 

• St • [T.(f (k ' 2) ) + T v (f m *) + D(f (0) ) - Q • D(f m * 

St • [T(f (k_1) ) + T v (f m *) + D(f <0) ) - Q • D(f m *) ] 
^m+1 _ ^(k) 


)] 


The Courant limit remains substantially unaltered. 

The dissipative terms are given as follows (ref. 2) 


and 
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D m = 1/(6 • St) • 5Z (f m - f m ) 

j = 1 3 




1/6 


z 

j=l 


m 

(p j 


m 

P 


) 


The subscript j refers now to every surrounding node involved in the 
finite volume approximation, and the superscript m refers to local time 
m • St. The terms referring to time m* • St are updated only at specific 
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Iterations and assumed constant between two updatings. The updating rate is 
taken equal to NV iterations, where NV is equal to 25, as a result of a 
numerical optimization. 

The dissipative term D is an approximation of second order differences. 
Without correction, the artificial viscosity coefficient is of the order of 
6r 2 / 6 1 , and the scheme is only first order accurate. When the viscosity is 
corrected and Q = 1 - 0(6r), the viscosity has a coefficient of the order 
of 6r 3 /$t, and the scheme is second order accurate. 

Cqi and Cq 2 are vectors of coefficients, i.e., the viscosity coefficients 
are dependent on this particular equation, for an improved accuracy/convergence 
ratio (ref. 2) . 

The four-stage Runge-Kutta algorithm adopted here replaces an Euler algo- 
rithm, i.e., a one-stage Runge-Kutta algorithm (refs. 1 to 3). Due to the 
enhanced stability, the total amount of the numerical viscosity is considerably 
reduced. This is the main advantage of the present time integration, since the 
total CPU time remains substantially unaltered, the reduction in the number of 
time steps due to convergence being balanced by the increase in the computa- 
tional time required by every time step. If we assign two work units to the 
evaluation of the residual (there are both viscous and inviscid contributions) 
and one work unit to the evaluation of the dissipative term and we define the 
efficiency of the scheme as the permitted CFL number divided by the number of 
work units, the one-stage scheme has an efficiency of about 2. 6/(5 + 2/NV), 
i.e. about the same efficiency. 

The time step is evaluated according to the classical CFL stability limit 
all over the computational domain. It is taken slightly smaller than the local 
CFL number in order to take into account the neglected stability limit due to 
the viscous terms. 

Due to the efficient pseudo unsteady solution, the method appears to be 
rather fast, while the explicit finite volume discretization allows ease of 
understanding and computer programming. 


RESULTS 

The turbulence model has been previously applied to the computation of 
backward facing step flows, in order to test the prediction capability of the 
low Reynolds number formulation when dealing with strongly recirculating flows 
(ref. 2). The computed length of the recirculation region generally compares 
fairly well with the experimental one. For a nearly incompressible flow, a 
Reynolds number (based on step height and inlet flow conditions) of 42 000 and 
an expansion ratio of 0.66, the predicted reattachment length is 7.50 against 
a measured value of 7.33 (ref. 2). Further calculations performed by modifying 
the Reynolds number, the inlet turbulence level, or the inlet length scale of 
the turbulence have shown the expected behavior. 

Sample calculations are performed here on a transonic turbine profile. 

The blade tested is a mean section of a gas rotor blade. The blade cascade is 
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shown in figure 2. The blade coordinates and experimental velocity distribu- 
tions are given in (ref. 5), while the experimental kinetic energy loss coeffi- 
cient distributions are given in (ref. 8). 

Calculations were preformed on a very fine computational grid, made up of 
141 pitchwise lines, 91 between leading and trailing edge, and 61 pseudo 
streamlines. About 20x10-* iterations are performed before convergence. The 
total CPU time on a VAX 8800 is less than 60 h. 

The inlet total temperature is Tqi ~ 280 K. The inlet total pressure var- 
ies more or less linearly with the outlet isentropic Mach number, the latter 
defined as fol lows 

M 2 . s = {2/ ( Y -1) • [P 2 /P 01 ) (1_y)/y -l.]} 1/2 

and poi ~ 1.5 bar at M 2 i s ~ 0.8, pp l ~ 2.75 bar at M 2 i s ~ 1.4. The inlet 
Reynolds number, based on inlet conditions and blade chord, carries linearly 
with the outlet isentropic Mach number in the same range, It is R e i ~ 300 000 
at M 2 i s ~ 1.4. Calculations were performed for an inlet flow angle and an 
outlet isentropic Mach number given as 

a] = 24° M 2 , i s = MO, 1.42 

In the comparison of predicted and measured / 5 / blade isentropic Mach 
number distributions at M 2 i s ~ 1-42 in figure 3, agreement appears to be very 
satisfactory. 

In a comparison of predicted and measured /8 / outlet kinetic energy loss 
coefficients, defined as 

c = [<p 01/ p 02 ><y-i>/y - u/n/2 • ( Y -1> • m* J s ] 

where Pq 2 is a mass averaged value over a blade pitch, the predicted values 
are (, ~ 0.047 at M 2 i s ~ 1 .10, C ~ 0.032 at M 2 i s ~ 1-42, while the measured 
values are C ~ 0.048 at M 2 } s ~ 1.10, C ~ 0.030 at M 2 i s ~ 1-42, thus result- 
ing in an experimental uncertainty in the predicted values very close to the 
experimental one. 

The good accuracy is due to the use of a relatively refined mathematical 
model and of a numerical integration with high spatial resolution and reduced 
numerical viscosity. 

If only the blade isentropic Mach number distributions has to be predic- 
ted, a simpler mathematical model and a less refined numerical integration can 
be used. If the two equation turbulence model is replaced by a mixing length 
model, the total enthalpy is assumed to be constant, the time integration is 
performed with a simple one step scheme, and if only a 65 by 19 computational 
grid is used, the isentropic Mach number distributions are predicted within an 
engineering accuracy, as represented in figures 4, 5, 9, and 10. In / 9 / the 
length scale of the turbulent motions is given by using a classical expres- 
sion, derived from the formulation proposed by Nikuradse for pipe flows and 
from the Van Driest damping factor, In / 1 0/ , it is simply taken proportional 
to the size of the spatial discretization. 



In the prediction of the outlet kinetic energy loss coefficient, these 
simplified solutions allow neither a qualitative nor a quantitative agreement. 
If a finer grid is used, the loss can be reduced, but the qualitative behavior 
cannot be properly reproduced. In calculations performed with a 79 by 39 
computational grid/1/, the predicted values are C ~ 0.054 at M 2 i s ~ 1-10, 

C, ~ 0.055 at M 2 i s ~ 1.42, while the measured values are C ~ 0.048 at 
M 2 i s ~ 1 • >0, C ~ 0.030 at M 2 is ~ 1.42. 


CONCLUSIONS 

The paper has presented a method for the solution of the Navier Stokes 
equations in transonic cascade flow fields. 

The work has been focused on the development of a low-Reynolds number and 
compressible K-e turbulence model. The turbulence model includes low Reynolds 
number terms, so that the equations are valid all over the laminar, transi- 
tion, and turbulent regions. Furthermore, the model includes a density gradi- 
ent term to better simulate variable density effects. 

The use of a four-stage Runge-Kutta algorithm allows one to significantly 
reduce the numerical viscosity due to the dissipative terms, thus leading to a 
better accuracy. 

First results of calculations compare favorably with experimental data. 
Even if uncertainties still remain, due to the limit of the numerical solution 
algorithm (influence of grid refinement, numerical viscosity, . .) the pro- 
posed method appears to be adequate for the study of steady transonic flows in 
turbine cascades. 
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